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1.0  INTRODUCTION 

As  aircraft  become  capable  of  sustained  flight  at  higher  altitudes,  their  rocket-propelled 
weapons  are  being  subjected  to  environments  at  launch  and  boost  which  were  not 
anticipated  when  these  weapons  were  designed.  Of  particular  importance  is  the  phenomenon 
of  plume-induced  separation.  Such  separation,  if  severe,  can  subject  control  surfaces  to  low- 
speed  recirculating  flows  which  may  seriously  degrade  the  effectiveness  of  the  surface.  A 
simple  calculation  may  determine  that  separation  could  not  occur,  for  instance,  if  the  nozzle 
exit  pressure  is  less  than  or  equal  to  free  stream.  However,  if  it  is  determined  that  the  nozzle 
is  highly  underexpanded,  i.e.,  the  jet-to-free-stream  static  pressure  ratio  is  greater  than  six, 
then  the  determination  of  whether  separation  is  occurring  to  such  an  extent  as  to  cause 
difficulties  requires  either  wind  tunnel  testing  with  a  real  plume  or  a  numerical  solution  of 
the  equations  that  closely  describe  the  flow.  This  report  details  the  development  of  two 
predictive  techniques:  one  a  component  method,  where  the  inviscid  and  viscous  parts  of  the 
flow  are  treated  separately  and  coupled  at  their  interfaces,  and  the  other  a  Navier-Stokes 
(N-S)  method. 

Figure  1  is  a  schematic  of  a  representative  afterbody  flow  with  plume-induced 
separation.  The  separation  occurs  because  the  afterbody  flow  must  compress  to  turn 
outward  around  the  plume.  When  the  adverse  pressure  gradient  is  larger  than  that  which 
enables  attached  flow,  the  afterbody  flow'  lifts  off.  When  the  detached  flow  reattaches  at  the 
plume,  it  turns  through  a  further  compression,  causing  a  secondary  shock  which  coalesces 
with  the  separation  shock  and  forms  a  lambda  shock  pattern.  The  plume  boundary  also 
turns  and  another  shock  is  formed  which  finally  coalesces  with  the  barrel  shock  of  the 
plume.  This  complex  of  shocks  is  known  in  component  methods  as  the  trailing  shock 
system.  Within  the  separated  layer,  a  low-speed  recirculating  flow  exists.  In  the  component 
method  this  region  is  usually  assumed  to  be  at  constant  pressure.  This  is  because  most 
measured  data  show  a  distinct  pressure  profile  shape  along  the  afterbody  that  increases 
abruptly  through  the  compression  and  reaches  a  level  that  is  fairly  flat  relative  to  the  initial 
compression.  The  pressure  changes  little  until  very  near  reattachment  on  the  jet  plume.  The 
mathematical  description  of  such  a  viscous  interaction  is  provided  by  the  Navier-Stokes 
equations. 

Although  computational  costs  are  declining,  routine  solution  of  the  Navier-Stokes 
equations  is  still  not  practical.  It  is  therefore  important  to  explore  the  capabilities  of  existing 
approximate  techniques,  determine  their  accuracy,  and  determine  their  extent  of 
applicability.  The  initial  phase  of  this  work  was  thus  directed  toward  developing  and 
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extending  an  existing  component  method  that  had  proved  to  work  well  in  predicting  base 
properties  where  separation  on  the  afterbody  was  not  present.  A  requirement  arose, 
however,  for  a  predictive  technique  for  vehicles  at  high  angles  of  attack.  The  component 
methods  are  not  applicable  in  these  flows,  because  they  require  an  external  criterion  for 
separation  that  is  based  on  singular  point  separation,  which  is  no  longer  valid  in  three- 
dimensional  flows  (Ref.  1).  This  led  to  the  modification  of  a  three-dimensional  N-S  solver 
specifically  to  solve  the  plume-induced  separation  problem.  This  work  with  the  N-S  solver 
allowed  it  to  be  brought  into  use  when  another  drawback  to  the  component  method  became 
apparent:  at  low  free-stream  Mach  number  and  high  jet-to-free-stream  pressure  ratios,  the 
external  flow  is  forced  to  turn  beyond  the  supersonic  turning  limit,  causing  a  detached 
normal  shock.  The  component  method  was  not  formulated  to  account  for  mixed  flows. 

This  report  is  of  two  parts.  The  first  part  presents  the  approach  to  developing  the 
component  method  into  a  tool  for  predicting  plume-induced  separation.  The  second  part  is 
devoted  to  the  modifications  required  of  a  three-dimensional  N-S  solver  to  make  it  useful 
for  predicting  plume-induced  separation  at  realistic  jet-to-free-stream  pressure  ratios. 

2.0  COMPONENT  METHOD 


2.1  GENERAL 

The  development  of  a  component  method  capable  of  predicting  plume-induced 
separation  is  based  on  previous  work,  Refs.  2-4,  which  had  generalized  and  extended  the 
basic  mixing  and  base  flow  theory  of  A.  J.  Chapman  and  H.  H.  Korst,  Refs.  5-6,  for 
predicting  base  properties.  This  generalized  theory  was  proved  capable  of  predicting  base 
flows,  taking  into  account  the  boundary  layer  existing  at  separation,  base  bleed,  and  the 
interaction  of  gas  streams  differing  in  chemistry  and  total  enthalpy.  Because  of  the  wealth  of 
literature  on  this  method,  this  report  will  dwell  only  on  the  changes  required  of  the  theory  as 
developed  in  Ref.  3. 

2.2  BASIC  COMPONENT  THEORY 

The  basic  component  method  based  on  Korst’s  work  contains  three  distinct  analyses:  an 
inviscid  analysis  is  used  to  determine  the  overall  flow  structure  and  the  pressure  field;  a 
second  analysis  generates  a  set  of  integral  conservation  equations  which  describe  the  viscous 
flow  into  and  out  of  the  base  region;  and  the  final  analysis  develops  the  mixing  theory  which 
produces  the  flow  property  profiles  across  the  shear  layers.  Each  of  these  analyses  may  be 
developed  independently.  They  are  coupled  in  an  iterative  process.  A  solution  to  the  base 
flow  problem  is  achieved  when  the  set  of  integral  conservation  equations  is  satisfied. 


6 


AEDC-TR-84-1  8 


At  any  step  in  the  iterative  process,  the  base  properties  serve  as  boundary  conditions  to 
the  inviscid  set  of  equations  and  to  the  mixing  model.  The  base  pressure  adjusts  the 
boundaries  of  the  converging  inviscid  streams,  the  mixing  theory  locates  the  mixing  profiles 
on  the  inviscid  boundaries,  and  the  integral  conservation  equations  balance  the  energies  and 
mass  rates  into  and  out  of  the  base  region.  The  residual  error  of  the  integral  conservation 
equations  drives  the  iteration  process. 

For  the  method  to  have  closure,  certain  key  streamlines  must  be  located  within  the  shear 
layers,  Figs.  2  and  3.  The  most  difficult,  historically,  to  determine  is  the  stagnation 
streamline,  which  determines  what  part  of  the  converging  shear  layers  will  turn  back  toward 
the  base  and  what  part  will  be  turned  downstream  to  experience  rccompression.  There  are 
many  theories  available  for  determining  this  streamline,  all  using  empirical  data  to  locale  the 
stagnation  point.  The  methods  of  Refs.  3  and  4  use  a  minimum  of  empirical  information  to 
develop  a  very  workable  theory  that  has  been  tesLed  without  change  over  a  wide  range  of 
flow  configurations,  including  flows  with  enthalpy  and  chemistry  differences  at  widely 
different  Mach  numbers. 

The  coordinate  systems  used  in  the  component  method  are  sometimes  contusing.  The 
general  coordinate  system  is  axially  symmetric  with  X  being  the  axial  coordinate  and  R  the 
radial  coordinate;  however,  the  mixing  layers  have  their  own  system.  The  coordinate  along 
the  longitudinal  direction  of  the  shear  layer  is  f;  l  =  0  at  the  point  of  separation.  Y  is  the 
coordinate  normal  to  f,  and  X  is  the  coordinate  along  the  slip  surface  where  the  two  streams 
converge  after  separation  (see  Fig.  3). 

2.3  EXTENSION  OF  THEORY 

The  analysis  of  Ref.  4  is  extended  in  a  straight  forward  manner  to  account  for  separation 
occurring  upstream  of  the  base-afterbody  juncture.  This  involves  a  modification  (o  the 
inviscid  flow  solver  (method  of  characteristics)  to  accommodate  the  separation  shock  wave, 
as  illustrated  in  Fig.  2,  and  ihe  inclusion  of  a  separation  criterion.  Hahn,  Ruppert,  and 
Mahal  (Ref.  7)  reviewed  existing  criteria  extensively  to  determine  when  separation  occurs. 
For  turbulent  free-interaction  separation,  the  authors  suggest  the  criterion  proposed  by 
Reshotko  and  Tucker 

M2/M|  -  0.762  (1) 

That  is,  separation  will  occur  when  the  ratio  of  the  Mach  number  downstream  of  the 
separation  shock  to  that  upstream  of  the  shock  is  equal  to  0.762.  This  criterion  is 
particularly  easy  to  apply;  however,  its  main  drawback  is  that  the  separation  on  a  boaltail 
afterbody  is  not  truly  of  the  free-interaction  type.  The  upstream  influence  ol  the  changing 
curvature  of  the  body  is  not  taken  into  account.  It  would  be  expected  that  separation 
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predictions  would  be  more  accurate  for  large  separations  where  the  point  of  separation  is 
well  up  on  the  body. 

2.4  MIXING  THEORY 

As  shown  in  Fig,  2,  the  separated  region  is  bounded  by  two  shear  layers  which  are 
assumed  to  converge  at  a  slipline.  The  slipline  is  determined  from  the  axially  symmetric, 
rotational  method  of  characteristics  (MOC).  The  method  of  characteristics  is  used  also  to 
determine  the  inviscid,  constant  pressure  boundaries  of  the  separated  region.  At  the  inviscid 
intersection,  point  I,  the  two  mixing  layers  are  assumed  to  have  known  velocity  profiles: 

*  =  y  F  1  +  erf(„  -  ,p)l  +  j  f,,/ne  "  1 ”  "  ^  df  (2) 

L  .1  N  0 


where  i',1/n  is  the  initial  boundary-layer  power-law  profile  and 


and 


<rY 
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yp  - 


od 
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where  a  is  the  mixing  parameter.  Equation  (2)  is  the  widely  used  error  function  profile 
distorted  by  the  initial  boundary  layer.  As  the  boundary-layer  thickness  goes  to  zero,  the 
second  term  becomes  negligible.  This  occurs  both  when  the  separating  boundary  layer  is  very 
thin  and  when  the  distance  from  separation  to  reattachment  is  large,  i.e.,  as  f  becomes  large. 
When  T?p  becomes  small,  the  mixing  layer  is  considered  fully  developed.  The  mixing 
parameter  a  is  determined  semiempirically.  The  method  of  Ref.  8  shows  it  to  be  well 
approximated  by 


a 


q  02di) 
GU 


(5) 


where  Ck  is  given  as  0.5085/<jo.  In  the  present  work  the  incompressible  <t0  is  taken  as  9.0. 


Each  profile  at  point  I  has  the  transverse  coordinate  77  normal  to  its  respective  inviscid 
boundary,  increasing  in  value  toward  the  high-speed  stream.  Following  Ref.  6,  the  Y  =  0 
point  is  located  by  a  momentum  balance  as 


Au 

l  gu2dY 


JYU  -  VM 

eu2dY  | 
yl  r 


(6) 


8 


AEDC-TH-84-18 


This  relation  is  solved  for  YM.  This  increment  is  the  distance  from  the  inviscid  boundary  to 
the  zero  point  of  the  profile.  Thus  this  relation  locates  the  profiles  relative  to  their  inviscid 
boundaries.  The  edges  of  the  mixing  zones  are  somewhat  arbitrary,  but  the  points  that  have 
worked  well  over  many  conditions  are  determined  by  the  following: 

1.  Yy  is  located  where  the  velocity  is  99.998895  percent  of  the  high-speed  value. 

2.  Yl  is  located  at  the  point  where  the  velocity  is  0.001105  percent  of  the  high¬ 
speed  value. 

Two  streamlines  must  now  be  identified  before  the  conservation  relations  for  the  base 
flow  may  be  presented.  The  first,  the  dividing  streamline,  is  defined  as 

j  L  GudY  |  =  j  LeudY  |  (7) 

°  l'=o  yd 

This  relation  is  solved  for  the  location  Yy  of  the  dividing  streamline.  It  is  that  location  in  the 
profile  outside  (“outside”  considered  as  toward  the  high-speed  edge)  of  which  the  mass  of 
the  just-separated  flow  is  accounted  for.  The  second  streamline  is  the  stagnating  streamline 
which  is  located  in  the  profile  at  Y$.  As  stated  earlier,  this  streamline  divides  the  mixing  layer 
into  the  parts  that  either  turn  back  toward  the  base  or  have  enough  energy  to  proceed 
downstream  through  recompression.  Locating  this  streamline  requires  a  theory  of 
recompression  as  outlined  in  the  following. 

2.5  RECOMPRESSION 

Korst’s  original  theory  assumed  that  the  stagnating  streamline  had  a  stagnation  pressure 
equal  to  the  static  pressure  downstream  of  the  reeompression  shock.  This  was  shortly  found 
to  be  incorrect  by  Korst  and  others.  (It  has  become  nearly  a  liturgical  requirement  of  any 
critical  discussion  of  Korst’s  theory  to  dismiss  his  theory  because  of  this  initial  assumption 
which,  of  course,  has  nothing  to  do  with  a  basic  theory  that  has  proved  its  usefulness  if  by 
nothing  more  than  its  longevity.)  The  actual  pressure  is  somewhat  less  than  the  peak 
pressure.  A  procedure  for  determining  this  pressure  and  thus  the  stagnating  streamline  was 
developed  in  Refs.  3  and  4  and  is  summarized  here. 

In  the  process  of  computing  the  inviscid  fields,  the  trailing  shock  system  is  also 
computed.  The  trailing  shock  system  shown  in  Fig.  2  also  determines  the  angle  of  the  slipline 
which  serves  as  the  convergence  surface  for  I  he  two  shear  layers.  Its  computation  is 
straightforward.  The  relations 


si  s: 
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and 


Psi  =  PS2 


(9) 


must  be  satisfied  along  the  slipline.  These  two  relations,  coupled  with  the  oblique  shock 
relations  for  the  two  streams,  arc  sufficient  for  determining  the  slipline  position  and  angle. 

The  recompression  is  assumed  to  stretch  from  X|  to  X2  (see  Fig.  3)  and  to  follow  the 
functional  form  (Ref,  3) 


p  “  p"  =  sin2  T (,/2)  (10) 

P,  -  PB  L'  '  X2  -  *1 .1 

where  the  variable  X  is  incremented  along  the  slipline.  The  location  X!  is  the  beginning  of 
recompression  and  X2  the  end.  (Note  that  recompression  is  assumed  to  occur  on  both  sides 
of  the  slipline,  each  recompression  process  independent  of  the  other  and  coupled  only  by  the 
common  static  pressure  of  the  base  region  and  the  common  pressure  downstream  of  the 
trailing  shock  system.) 

The  location  of  X|  is  found  geometrically  by  making  the  angles  (3\  and  equal.  The 
location  of  X2  is  found  with  the  expression 

P  x2  _  f  v  u  +  Ym 

I  pRdX  -  sin7R  l  gu-dY  +  pB(Ru  +  R,)  (Xu  -  X,)/2  (11) 

x,  y, 


which  is  developed  in  Ref.  3  from  a  momentum  balance  of  the  oncoming  shear  layer  with  the 
pressure  force  on  the  slipline.  The  radial  effect  has  been  taken  into  account  because  of  the 
sometimes  significant  change  in  radius  along  the  slipline  (the  radial  effect  in  the  mixing  layer 
is  relatively  minor).  The  stagnating  streamline  is  assumed  to  come  to  rest  isentropically. 
Thus  sufficient  information  is  developed  to  determine  the  location  of  the  stagnating 
streamline.  For  example,  an  isoenergetic  perfect  gas  gives 


1  ~  (iW’pJb  1J/~' 

1  -  (Pb/Ptu)1-1  "  1);> 


02) 


The  geometry  of  the  profile  relative  to  the  slipline  and  Eqs.  (10)  -  (12)  completely  determine 
the  location  Ys  of  the  stagnating  streamline. 


2.6  CONSERVATION  EQUATIONS  OF  THE  BASE  FLOW 


With  the  key  streamline  locations,  Y\|,  Yd.  and  Ys  known,  it  is  possible  to  set  down  the 
basic  mass  and  energy  balance  relations.  For  the  conservation  of  mass, 


10 


AEDC-TR  84-18 


VL 

and  for  energy, 
Ys 


fY° 

\  eudY 


i  Yd 


YS 


I  +  f  pudY  I  -  j  pudY 
si  yl  s2  yl 


Sl 


rVs 

\  pudY 

yl 


(13) 


S2 


f  euHdY  | 
S) 


+  f  puHdY  I  =  Hb  (  (  e^dY  +  (  e^Y  I  ^ 
i  L  VC  <n  *>• 

These  constitute  the  primary  equations  which  are  solved  for  the  properties  in  the 
base/separated  region.  Several  auxiliary  relations  may  be  used  when  a  species  difference 
exists  between  the  exhaust  gas  and  the  free  stream.  These  are  found  in  Ref.  3  and  are  not 
detailed  here. 


2.7  INITIAL  VELOCITY  PROFILE 


The  theory  as  developed  in  Refs.  3-4  assumed  that  the  separating  layer  off  a  bluff  base 
would  have  a  power-law  profile.  For  separation  occurring  upstream  of  the  base,  the 
separating  profile  is  more  closely  related  to  the  wake-like  profile  of  a  developing  shear  layer 
as  shown  in  Fig.  4,  using  data  from  Refs.  9  -  10.  This  shape  develops  in  a  shear  layer  from  an 
initial  power-law  profile  some  distance  downstream  of  a  bluff  edge.  To  more  closely  model 
the  shape  of  the  separating  profile,  the  separating  layer  is  assumed  to  separate  with  a  profile 
with  tjp  =  4.0.  Thus,  in  the  present  calculations,  it  is  assumed  that  the  distance  f  has  been 
extended  by  a  virtual  amount.  Thus  the  effective  mixing  length  is 


4ff  =  f  +  fnr 

where 

4 =  05/4.0  06) 

2.8  COMPUTATIONAL  PROCEDURE 

Initial  calculations  were  done  using  the  GASL  method  of  characteristics  program,  Ref. 

1 1.  This  code  was  found  to  convect  total  pressure  loss  from  the  bow  shock  too  far  into  the 
field  of  the  afterbody.  This  caused  an  overprediction  of  the  extent  of  separation  at  higher 
Mach  numbers.  The  present  results  were  calculated  by  assuming  a  constant  pressure  equal  to 
the  free-stream  pressure  from  the  shoulder  of  the  missile  to  just  upstream  of  the  boattail 
break.  Boundary-layer  growth  was  then  determined  over  this  length  using  the  method  of 
Ref.  12.  The  solution  procedure  is  then  as  follows  for  isoenergetic  flows: 
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1.  The  MOC  of  Ref.  3  is  started  assuming  free-stream  conditions  with  no 
boundary-layer  displacement  taken  into  account.  This  computation  is  marched 
downstream  to  an  assumed  point  of  separation. 

2.  A  separation  shock  wave  is  computed  and  the  MOC  is  continued.  The 
boundary  on  the  edge  of  the  separated  region  is  determined  by  the  guess  of 
pressure  in  that  region.  Computation  of  the  plume  field  is  also  done  with  the 
MOC.  The  intersection  of  the  plume  and  outer  stream  is  found  and  used  to 
determine  the  length  of  the  inviscid  boundaries. 

3.  The  trailing  shock  system  and  slip  surface  are  computed  to  obtain  the  peak 
recompression  pressure. 

4.  The  viscous  mixing  theory  is  then  applied  and  the  key  streamlines  are 
determined.  Finally,  a  mass  balance  of  flow  into  and  out  of  the  separated 
region  is  made. 

If  the  mass  balance  equation  is  not  satisfied,  an  error  is  determined,  a  new  pressure  is 
calculated,  and  steps  2  through  4  are  repeated.  Once  a  converged  separation  pressure  is 
determined,  the  separation  criterion  for  a  turbulent  boundary  layer,  Eq.  (1),  is  used  to 
determine  whether  separation  has  occurred.  If  this  condition  is  not  satisfied,  a  new 
separation  point  is  assumed,  and  steps  1  through  4  are  repeated.  This  continues  until  a 
separation  point  is  determined. 

2.9  DISCUSSION  OF  COMPONENT  METHOD  RESULTS 

Experimental  data  for  afterbodies  with  a  plume-induced  separation  are  not  widely 
available.  However,  two  reports  of  measured  data  that  have  appeared  in  the  literature  are 
sufficient  for  validating  the  method  as  an  engineering  tool  (Refs.  13  and  14).  The  results 
from  Ref.  13  are  from  a  configuration  run  at  a  free-stream  Mach  number  of  3.5  at  various 
static  jet-to-free-stream  pressure  ratios.  Comparisons  are  presented  in  Fig.  5  of  predictions 
made  using  the  present  component  method  with  the  measured  data.  Actual  measured 
separation  points  are  not  presented  as  they  were  not  available;  however,  the  predicted 
separation  point,  which  is  indicated  by  the  step  rise  in  pressure,  is  located  at  a  position 
within  the  pressure  rise  that  is  consistent  with  the  beginning  of  separation  from  other  data 
that  are  available,  i.e..  Ref.  14.  They  show  separation  occurring  10  to  20  percent  into  the 
pressure  rise.  The  calculated  plateau  pressure  is  high  compared  to  the  measured  data.  This  is 
because  the  predicted  pressure  is  that  which  occurs  downstream  of  an  inviscid  oblique  shock 
wave,  which  is  of  zero  thickness.  The  measured  pressure  on  the  afterbody,  however,  is  the 
result  of  the  smearing  effect  of  the  viscous  layer,  which  spreads  out  the  compression.  The 
peak  pressure  is  thus  pushed  downstream,  perhaps  off  the  body. 
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Comparison  was  made  with  the  data  of  Ref.  14  at  a  lower  Mach  number.  This 
comparison  is  shown  in  Fig.  6.  At  the  lower  pressure  ratios,  the  comparison  is  quite  poor.  In 
general,  the  component  method  gives  quite  acceptable  results  at  high  jet-to-free-stream 
pressure  ratios  where  the  separation  is  greater  than  three  nozzle  radii  up  the  afterbody  from 
the  base.  This  is  expected  since  the  external  criterion  used  to  determine  when  separation  is 
occurring  is  based  on  free  interaction  separation  where  separation  is  based  only  on  local 
conditions  free  from  direct  influences  of  downstream  geometry  such  as  the  afterbody/base 
juncture. 

An  anomalous  condition  sometimes  occurs  when  trying  to  predict  separation  on  bodies 
where  the  Mach  number  of  the  outer  stream  is  significantly  lower  than  the  Mach  number  of 
the  expanded  jet  boundary  at  the  slipline.  If  the  pressure  ratio  of  the  jet  to  free  stream  is 
sufficiently  high,  then  a  solution  to  the  trailing  shock  problem  cannot  be  obtained  because 
the  outer  stream  is  required  to  turn  further  than  the  oblique  shock  limit.  That  is,  the 
recompression  shock  on  the  plume  boundary  must  be  detached.  This  means  that  a  subsonic 
region  exists  downstream  of  the  detached  shock  wave.  The  basic  theory  of  this  component 
method  does  not  account  for  such  flows;  therefore,  this  type  of  flow  must  be  treated  by 
other  methods,  i.e.,  Navier-Stokes  methods.  This  condition  played  substantially  in  the 
decision  to  pursue  the  Navier-Stokes  methods  as  discussed  in  the  following. 

3.0  THE  NAVIER-STOKES  APPROACH 


3.1  GENERAL 

There  are  two  principal  reasons  for  exploring  the  use  of  N-S  solvers  for  the  solution  of 
the  plume-induced  separation  problem: 

1.  The  occurrence  of  the  detached  shock  at  the  plume  boundary  with  its  attendant 
subsonic  region  invalidates  the  component  method  because  of  its  use  of  a 
spatially  hyperbolic  inviscid  solver  (MOC).  This  happens  when  the  two  streams 
are  of  substantially  different  Mach  numbers. 

2.  A  requirement  for  predictions  at  very  high  angles  of  attack  (10  to  20  deg)  also 
invalidates  the  component  method,  because  its  externally  applied  separation 
criterion  is  based  on  the  occurrence  of  singular  type  separation,  which  does  not 
occur  in  general  in  three-dimensional  flows. 

A  third,  but  less  justifiable,  reason  is  the  fact  that  N-S  solvers  are  conceptually  pleasing 
because  they  use  fewer  empirically  determined  constructs.  That  element  of  the  method 
which  does  require  empiricism,  the  modeling  of  the  turbulent  viscosity,  shows  signs  of 
becoming  a  rapidly  maturing  technology  (Refs.  15  -  16). 
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When  the  N-S  solver  was  used  as  detailed  here,  it  became  apparent  early  that  extensive 
modification  of  the  code  would  be  required,  if  the  extent  of  separation  at  any  jet-to-free- 
stream  pressure  ratio  of  interest  was  to  be  computed.  The  difficulty  arises  in  the  region  just 
downstream  of  the  nozzle  lip.  At  high-pressure  ratios,  the  barrel  shock  of  the  plume  begins 
to  form,  developing  a  region  of  very  high  gradients.  The  expansion  from  the  nozzle 
overshoots,  causing  the  local  pressure  to  be  lower  than  the  ambient  pressure.  The  flow  must 
then  be  compressed  sharply  through  the  mechanism  of  the  barrel  shock  to  the  boundary 
pressure.  The  usual  method  of  overcoming  this  type  of  problem  is  to  either  concentrate  the 
computational  mesh  in  this  region  or  add  additional  smoothing  or  do  both.  The  adding  of 
points  quickly  reaches  a  limit  when  the  available  memory  in  the  computer  is  exhausted  and 
the  resolution  of  other  important  features  of  the  flow  is  seriously  degraded.  Smoothing  fails 
for  an  interesting  reason.  Smoothing  effectively  increases  the  viscosity  in  the  high  gradient 
region.  This  has  the  effect  of  over-entraining  fluid  from  the  region  of  the  base  just  above  the 
nozzle  lip  to  the  point  that  the  pressure  becomes  unrealistically  low,  often  to  the  point  of 
being  negative,  causing  the  time-marching  scheme  to  fail.  If  heroic  measures  are  taken  to 
prevent  the  low  pressure,  the  converged  solution  will  often  lack  sufficient  predictive 
accuracy.  For  these  reasons,  an  artifice  was  developed  to  circumvent  the  nozzle-lip  problem. 

Because  of  the  behavior  of  the  algorithm  in  the  nozzle-lip  region  and  from  experience 
derived  from  development  of  the  component  model  discussed  earlier,  the  interaction  at  the 
nozzle  lip  was  determined  to  be  primarily  inviscid.  The  nozzle  lip  is  thus  treated  inviscidly 
and  removed  from  the  region  of  the  viscous  solver.  The  inviscid  MOC  is  used  to  develop  new 
boundary  conditions  for  the  N-S  solver  downstream  of  the  nozzle  lip. 

Results  are  presented  of  comparisons  with  available  experimental  data  using  this  new 
procedure.  A  study  using  the  same  measured  data  for  comparison  was  done  by  Deiwert, 
Ref.  17.  He  used  an  azimuthally  invariant  form  of  the  equations  and  a  different  gridding 
philosophy. 

In  the  following,  since  the  numerical  tools  are  well  developed  and  readily  available  in  the 
open  literature,  many  of  the  details  will  be  referenced. 

3.2  BASIC  EQUATIONS 

The  three-dimensional,  unsteady,  thin-layer  Navier-Stokes  equations  are  transformed 
(see  Appendix  A)  into  the  curvilinear  coordinate  system  £,  77,  f,  and  7  as 

dA  +  d£  +  3,F  +  drG  =  Re'>  3fS  (17) 
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where,  with  reference  to  Fig.  1,  £  is  the  body-conforming  coordinate  in  the  streamwise 
direction,  y  is  the  body-conforming  coordinate  in  the  azimuthal  direction,  and  f  is  the 
outward  coordinate  ray  from  the  body  to  the  outermost  boundary.  The  Cartesian  x 
coordinate  is  the  axis  of  the  body,  the  z  coordinate  is  normal  to  x  and  tangent  to  the  base  of 
the  body,  and  the  x-z  plane  forms  the  slice  used  for  presentation.  Thus  the  Cartesian  velocity 
components  of  interest  are  u  and  w.  The  vectors  q,  E,  F,  and  G  are 

GU 

GiiU  +  £xP 

E  =  J-'  gvU  +  £jp  U8> 

qwU  +  £ZP 
(e  +  p)U  -  £tP 

qV  eW 

guV  +  r?xP  „  GuW  +  f*P 

F  =  J-l  gvV  +  rjyp  G  =  J_i  GvW  +  fyp 

gwV  +  i?zp  GwW  +  fzp 

(e  +  p)V  —  ijtP  (e  +  P)W  -  ftP 

0 

+  fy  +  fl)  Uf  +  (^/3Xr,Uf  +  fyVJ  +  fzWf)fx 

+  rD  vf  +  (/t/3)Cfxur  +  ryvf  + 

S  =  J-l  +  f2  +  f2)  Wf  +  +  fyvr  + 

[(r5  +  Ty  +  f|)[o.5^(u2  +  v2  +  w2)f 

+  kPr-i(7-  1)~!  (a2)r]  +  (|*/3Xf*u  +  fyv  +  f4w) 

X(fxuf  +  fv  vr  +  ^7wf)] 

and  the  contravariant  velocities  are  _ 

U  =  it  +  +  izw 

V  =  rjt  +  rjxU  +  i)yv  +  rjzw  (21) 

w  =  r,  +  rxu  +  fyv  +  rzw 
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The  pressure  is  determined  from 

p  =  (7  -  l)[e  -  0.5g(u2  +  v2  +  w2)]  (22) 

and  the  Cartesian  velocity  components  are  non-dimensionalized  with  the  free-stream  speed 
of  sound  a^;  density  q  is  scaled  by  and  the  total  energy,  e,  by 

The  chain  rule  expansion  of  derivatives  of  the  Cartesian  coordinates  with  respect  to  the 
curvilinear  coordinates  is  solved  to  give  the  metric  terms 


£x  =  J(y,zr  -  y^v) 

£y  =  J(z^f  -  X,z f) 

it  =  J(x„yt  -  y^r) 
f*  =  J(y^  -  W 

fy  =  -  XfZ^) 

tz  =  J(x£y„  -  y**,) 


Vx  =  Jtejyr  -  yjZf) 

T?y  =  -  xfZ£) 

Vx  =  Hypt  -  x£yr) 

it  =  -  *r£x  -  y^y  “  Zrfz 
*Jt  =  ^tVx  yr^y  —  ZT)Jz 
ft  =  —  Xrfx  —  yrfy  ~  ZTfz 


(23) 


with 

J-i  =  x*  y,  zr  +  xf  y£  z,  +  x,  y s  z£  -  x£  yr  z,  -  x,  y£  zf  -  xf  y,  z£  (24) 
The  approximate  factorization  difference  equation  is 

(1  +  h6£  An  -  e,  J“L  V£A£J)(I  +  ha,  Bn  -  ei J " 1  V,A,J)x 
(I  +  hfif  Cn  -  hRe-‘  5f  J~ 1  MnJ  -  ei  J"1  VrArJ)(qn+l  -  qn) 

(25) 

=  -  At  (af  En  +  5,  Fn  +  af  Gn  -  Re-*  6fSn) 

-  cEJ->  [(V£A£)2  +  (V,A,)2  +  (VrAr)2]  Jq" 

where  the  <5’s  are  the  central-difference  operators,  A  and  V  are  forward  and  backward 
difference  operators,  and  h  =  t/2.  The  matrices  An,  B",  and  Cn  are  obtained  from  the 
linearization  in  time  of  En,  Fn,  and  Gn.  These  are  detailed  by  Pulliam  and  Steger,  Ref,  18, 
along  with  the  coefficient  matrix  Mn.  The  smoothing  terms  of  the  forms  e|  J  - 1  V£A£  Jq  and 
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eE  J~'  (V^)*  Jq  are  added  specifically  to  damp  nonlinear  instabilities  of  the  central 
difference  scheme. 

This  algorithm  has  been  put  into  a  practical  code  by  Pulliam  and  Steger.  The  actual  code 
used  herein  is  a  vectorized  form  of  Pulliam  and  Steger’ s  code  developed  by  J.  A.  Benek. 

3.3  GRID 

The  development  of  a  grid  for  the  computational  domain  is  perhaps  as  important  as  the 
solution  algorithm  itself.  Because  the  thin-layer  form  of  the  conservation  equations  is  being 
used,  certain  physical  assumptions  are  carried  with  any  choice  of  grid. 

The  principal  assumption  made  in  the  choice  of  grid  is  that  the  wake  flow  downstream  of 
the  base  region  including  the  plume  may  be  treated  as  an  inviscid,  rotational  flow.  It  is  also 
assumed  that  the  boundary  layer  in  the  nozzle  has  little  effect  on  separation  and  may  be 
neglected.  The  effluent  from  the  nozzle  is  thus  treated  as  an  inviscid,  conical  flow. 

These  basic  assumptions  lead  to  the  backward  “c”  grid  as  shown  in  Fig.  7.  By  wrapping 
the  afterbody/base  region  in  a  body-conforming  coordinate,  the  basic  code  is  left 
unchanged;  concentrations  of  rays  may  be  easily  positioned  at  the  nozzle  lip  region  where 
resolution  is  most  important,  and  the  rays  can  be  readily  aligned  with  the  plume  boundary  so 
that  the  reattachment  point  and  the  trailing  shock  system  may  be  better  resolved.  Since  the 
thin-layer  assumption  retains  only  those  viscous  terms  in  the  conservation  equations  that 
have  derivatives  in  the  f  direction  and  since  the  f  coordinate  is  nearly  aligned  with  the 
principal  flow  direction,  the  viscous  terms  effectively  eliminated  are  in  the  region 
downstream  of  the  base/afterbody  juncture. 

The  grid  is  wrapped  in  opposite  directions  away  from  the  nozzle  lip  with  an  exponentially 
increased  spacing.  A  radius  equal  to  10  percent  of  the  nozzle  radius  is  added  at  the  juncture 
of  the  base  and  the  afterbody  to  smooth  the  grid  around  this  corner.  Axially  symmetric  flow 
is  enforced  by  using  five  planes  with  identical  boundary  conditions.  Five  planes  are  required 
because  of  the  fourth-order  smoothing  required  for  stability  of  the  algorithm,  Eq.  (25). 
Early  experimentation  with  underexpanded  nozzles  led  to  a  device  to  avoid  the  singularity  at 
the  nozzle  lip.  Without  extremely  dense  packing  of  points  in  the  lip  region,  no  flow  with 
pressure  ratios  of  significance  could  be  computed.  Consequently,  the  nozzle  lip  was 
eliminated  from  the  region  of  computation  for  the  viscous  code.  This  was  accomplished  by 
adding  a  new  surface  extending  from  just  above  the  nozzle  lip  to  the  nozzle  centerline,  as 
shown  in  Fig.  8. 
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3.4  BOUNDARY  CONDITIONS  AND  NOZZLE-LIP  ARTIFICE 

The  Pulliam-Steger-Benek  N-S  solver  uses  explicit  boundary  conditions.  While  it  is 
arguable  that  this  type  boundary  condition  impedes  convergence,  it  nevertheless  provides  an 
extremely  practical  code  in  which  boundary  conditions  may  be  easily  input  or  modified  as 
were  done  on  the  new  surface  that  avoids  the  nozzle-lip  problem. 

The  boundary  conditions  on  the  new  surface  are  provided  by  the  MOCs  (see  Ref.  3),  as 
shown  in  Fig.  9.  Beginning  with  the  exit  plane  conditions,  the  steady,  Euler-equations  MOC 
is  used  to  march  an  inviscid  field  out  to  the  new  boundary.  The  nozzle-lip  singularity  is 
treated  as  a  multivalued  point  using  a  Prandtl-Meyer  centered  expansion.  The  flow  at  the  lip 
is  expanded  until  it  reaches  the  pressure  of  the  base  point  that  is  the  intersection  of  the  new 
boundary  and  the  base.  As  the  base  pressure  changes  continuously  from  initial  conditions  to 
the  final  converged  solution,  the  expansion  is  adjusted  with  each  time  step.  The  boundary 
conditions  of  the  Navier-Stokes  solver  are  thus  being  adjusted  with  each  step.  The  boundary 
conditions  on  the  new  surface  that  are  imposed  on  the  Navier-Stokes  solver  are  obtained 
from  the  characteristics  by  a  double  linear  interpolation:  i.e.,  the  first  interpolation  is 
performed  as  the  characteristics  cross  the  new  boundary,  and  the  second  is  done  to  find  the 
flow  properties  at  the  fixed  grid  of  the  Navier-Stokes  solver.  The  MOC  develops  its  own  grid 
from  the  initial  spacing  in  the  nozzle  exit  plane  and  from  the  number  of  pressure  decrements 
through  the  expansion.  So,  a  finer  grid  spacing  in  the  nozzle  exit  plane  or  in  the  expansion 
will  increase  the  resolution  on  the  boundary  for  the  N-S  solver.  A  finer  expansion  grid  is 
required  for  the  higher  pressure  ratios. 

The  boundary  conditions  elsewhere  are  imposed  in  a  straightforward  manner.  On  the 
body,  the  no-slip  conditions  are  enforced  with  the  contravariant  velocities  set  to  zero. 
Density  is  imposed  on  the  body  as  a  first-order  extrapolation  from  the  first  grid  point  off  the 
body.  Adiabatic  wall  conditions  are  enforced  with  the  total  energy,  e,  determined  from  the 
zero  pressure  gradient  condition  normal  to  the  wall.  The  upstream  boundary-layer  profile  is 
taken  from  the  experimental  data.  Free-stream  conditions  are  imposed  at  the  outer 
boundary  from  the  upstream  boundary  to  the  coordinate  ray  through  the  afterbody/base 
juncture,  where  the  condition  is  changed  to  a  zero-gradient  outflow  condition  from  that 
point  to  the  centerline.  The  conditions  between  the  last  streamline  of  the  expansion  fan  and 
the  base  are  obtained  by  first-order  extrapolation  from  the  field. 

3.5  TURBULENCE  MODEL 

For  N-S  solvers  to  be  useful  in  predicting  flow  behavior  for  practical  configurations,  they 
require  a  turbulence  model  to  relate  properties  of  the  flow  field  to  the  apparent  change  in 
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in  viscosity  caused  by  turbulence.  The  algebraic  models  are  most  attractive  because  of  their 
simplicity. 

Baldwin  and  Lomax,  Ref.  15,  developed  an  algebraic  turbulent  viscosity  model  that  is 
widely  used  because  of  its  favorable  trade-off  between  ease  of  use  and  predictive  accuracy  in 
unseparated  flows.  This  model  is  a  two-layer  model  with  the  inner  mixing  length 
proportional  to  the  product  of  distance  from  the  wall  and  the  Van  Driest  damping  factor.  It 
crosses  over  to  a  wake-type  model  in  which  the  mixing  length  is  scaled  by  the  vorticity. 
According  to  Ref.  15,  the  model  is  valid  for  separated  flows  and,  with  the  outer  model 
alone,  for  pure  wakes.  However,  Thomas  (Ref.  19)  reports  that  if  the  outer  model  alone  is 
used  for  wakes,  instabilities  develop,  and  that  Baldwin,  in  private  communication  with 
Thomas,  recommended  using  the  basic  two-layer  model  with  minor  adjustment.  Thomas 
proceeded,  however,  to  develop  his  own  variant.  But  for  the  present  work,  it  was  determined 
from  numerical  experiments  that  the  basic  Baldwin-Lomax  two-layer  model,  with  minor 
change,  worked  quite  well. 

The  eddy  viscosity  is  defined  as 


£in  =  e*2  I  “i  (26) 

near  the  body.  It  is  defined  in  the  wake  as 

eout  =  KCcp  g  Fwake  Fkleb  (y)  (27) 

The  inner  model  is  switched  over  to  the  wake  model  when  qn  =  rout.  The  variable  y  here 
refers  to  normal  distance  from  the  body  and 


and 


®  i  =  [(dyu  -  dxv)2  +  (azv  -  3yw)2  +  (axw  -  dzu)2] 
f  -  ky[l  -  exp  (-J71)] 


1/4 


where 


Also, 


y+  =  y 


[ 


wall 


(28) 


(29) 


(30) 


(31) 
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and 


F 


wake 


y max  F, 


max 


(32) 


or 


F«,ake  —  Ymax  U^jf  Fn 


(33) 


depending  on  which  is  the  smaller.  The  terms  yma*  and  Fmax  are  determined  from  the 
maximum  of  the  function 


Also, 


F(y)  =  y 


w 


[,  „p(^L)] 


Udif  =  (u2  +  v2  +  w2)l/j  -  (u2  +  v2  +  w2)/2. 


(34) 


(35) 


Baldwin  and  Lomax  assigned  the  following  values  to  the  parameters: 


A+  =  26.0 

CCp  =  1.6 

Ckleb  =  0*3 
Cwk  =  0.25 
k  =  0.4 
K  =  0.0168 


(36) 


3.6  TURBULENCE  MODEL  VERIFICATION 

Numerical  experiments  were  performed  for  the  configuration  shown  in  Fig.  10  to 
evaluate  the  two-layer  model  in  separated  flow.  This  configuration  of  an  afterbody  and  solid 
plume  simulator  was  tested  in  the  Acoustic  Research  Tunnel  (ART)  at  AEDC,  Ref.  20.  As 
shown  in  Fig.  10,  the  pressure  distribution  and  the  separation  point  prediction  were 
improved  significantly  if  it  were  assumed  that  the  boundary  layer,  after  separation,  became 
essentially  free  of  turbulent  stress.  A  similar  observation  was  made  by  Swanson,  Ref.  21, 
relative  to  solutions  of  the  Navier-Stokes  equations  about  a  configuration  very  similar  to  the 
ART  model.  He  found  that  a  relaxation  turbulence  model  suggested  by  Shang  and  Hankey, 
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Ref.  22,  improved  his  results  significantly.  According  to  Swanson,  this  was  because  the  eddy 
viscosity  was  decreased  in  the  outer  (separated)  layer,  which  served  to  decrease  the  plateau 
pressure,  bringing  it  into  better  agreement  with  measurements.  Good  results  were  obtained 
with  the  ART  model  by  restricting  the  turbulent  viscosity  to  the  region  near  the  body  no 
larger  than  the  attached,  upstream,  in-flow  boundary-layer  thickness.  As  a  practical  matter, 
this  was  accomplished  by  terminating  the  calculation  of  turbulent  viscosity  at  the  grid  point 
away  from  the  body  surface  nearest  the  thickness  of  the  boundary  layer.  A  family  of 
pressure  distribution  curves  between  the  upper  and  lower  curves  of  Fig.  10  may  be  calculated 
by  varying  the  point  of  termination  of  turbulence  from  the  thickness  of  the  boundary  layer 
to  the  width  of  the  entire  field.  No  calculations  were  made  with  the  turbulence  turned  off 
any  closer  to  the  body  than  the  boundary-layer  thickness.  Because  of  these  Findings,  the 
turbulence  was  calculated  only  out  to  the  thickness  of  the  upstream  boundary  layer  in  the 
comparisons  with  the  FFA  measurements  that  follow. 

Since  these  numerical  experiments  with  the  Baldwin-Lomax  model  were  conducted,  a 
new  model  has  been  introduced  by  Johnson  and  King,  Ref.  16,  that  promises  to  overcome 
the  above  difficulties.  The  model  predicts  more  accurately  the  apparent  viscosity  in  the  outer 
region  of  the  separating  boundary  layer  and,  according  to  the  paper,  gives  excellent 
predictions  of  the  pressure  distribution  in  the  separated  region. 

3.7  COMPUTED  RESULTS  AND  COMPARISONS  WITH  MEASURED  DATA 

The  method  was  applied  to  a  configuration  run  by  Agrell  and  White,  Ret.  14,  as  shown 
in  Fig.  11,  where  the  unit  Reynolds  number  is  I95,000/RN,  Rn  is  15  mm,  and  the  afterbody 
boundary-layer  shape  is  the  1/7  power-law  profile. 


The  configuration  with  a  boattail  of  8  deg  was  chosen  because  it  was  the  one  for  which 
the  most  results  were  shown  in  Ref.  14.  The  initial  conditions  used  were  free  stream 
everywhere  with  the  upstream  boundary-layer  profile  and  thickness  assumed  along  the  entire 
body-conforming  coordinate.  The  nozzle  exit  conditions  were  set  to  the  correct  Mach 
number  for  a  conical  How  nozzle,  but  the  exit  static  pressure  was  initially  set  equal  to  the 
ambient  pressure.  After  approximaely  300  time  steps,  the  nozzle  exit  pressure  was  increased 
in  stages  up  to  the  static  jet-to-free-stream  pressure  ratio  of  interest.  The  maximum  pressure 
ratio  of  the  experiments  was  1 5.2.  No  upper  limit  to  pressure  ratio  was  determined;  however, 
the  higher  pressure  ratios  were  more  difficult  to  obtain  because  smaller  time  steps  had  to  be 
used  to  assure  convergence.  Achieving  a  higher  pressure  ratio  using  a  converged  solution  as  a 
starting  point  was  usually  more  difficult  than  beginning  from  the  initial  conditions  and 
stepping  the  pressure  every  few  hundred  time  steps. 
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The  results  and  comparisons  are  presented  in  Figs.  11-20.  The  predicted  point  of 
separation  on  the  afterbody  is  compared  with  the  measured  results  in  Fig.  11.  The 
comparison  is  very  good  at  the  lower  pressure  ratios  with  a  gradual  underprediction 
noticeable  as  the  pressure  ratio  passes  through  twelve.  This  underprediction  is  attributed  to 
the  decrease  in  grid  resolution  on  the  plume  boundary  at  the  reattachment  point,  because  an 
increase  in  the  pressure  ratio  lifts  the  plume  boundary  out  of  the  finer  mesh  region.  This 
situation  could  be  helped  by  using  an  adaptive  grid. 

Also  in  Fig.  1 1  is  the  separation  curve  computed  using  the  Chapman-Korst  component 
method  as  described  previously  and  shown  in  Fig.  6.  It  shows  poor  agreement  for  the  lower 
pressures.  This  is  caused  by  the  method’s  reliance  on  theory  developed  for  free-interaction 
separation  before  a  forward  facing  step.  This  theory  does  not  account  for  the  effect  of  the 
afterbody  shape  with  bluff  base  on  separation.  Thus  when  the  separation  becomes  quite 
large,  the  base  effects  being  propagated  upstream  are  diminished  and  agreement  improves. 
The  Chapman-Korst  component  method  is  used  at  the  AEDC  as  a  quick  screening  method 
to  detemine  whether  separation  could  be  a  problem. 

Figure  12  presents  the  comparison  with  experiment  of  the  afterbody  pressure 
distribution.  The  predicted  pressure  distribution  shows  good  agreement,  except  for  the 
plateau  region.  However,  if  the  first  peak  is  compared  with  the  experimental  plateau 
pressures  as  in  Fig.  13,  the  agreement  is  good.  The  dip  in  the  pressure  curve  after  the  first 
peak  appears  to  be  attributable  to  the  lack  of  grid  resolution  and  to  the  radius  inserted  in  the 
base/afterbody  juncture. 

The  base  pressure  shows  poor  agreement  with  measured  data  as  shown  in  Fig.  14,  The 
base  pressure  has  a  fairly  large  variation  on  the  bluff  face.  This  is  caused  both  by  the  base 
corner  radius  and  by  the  displaced  boundary  at  the  nozzle  lip.  No  experimentation  was 
performed  to  determine  the  smallest  radius  allowable  by  the  method.  It  is  apparent  from 
calculations  made  on  different  configurations  that  this  radius  is  much  too  large  for  the  base 
size.  This  size  has  caused  the  recirculating  eddies  to  be  displaced,  causing,  in  turn,  the 
pressure  disparity.  Calculations  made  with  different  configurations  also  show  that  if 
sufficient  flat  base  area  is  maintained  relative  to  the  corner  radius  and  the  displaced 
boundary,  a  typical  flat  pressure  distribution  on  the  face  is  obtained.  Experimentation  with 
the  position  of  the  intersection  of  the  displaced  boundary  and  the  base  area  showed  that 
moving  the  point  changed  the  local  base  pressure,  but  not  the  position  of  separation. 

The  contour  plots  of  density  and  static  pressure  in  Figs.  15-18  show  the  computed 
features  of  the  entire  field.  The  lambda  shock  structure  is  readily  apparent,  and  it  is 
particularly  sharply  defined  in  the  pressure  contour  plot  at  a  pressure  ratio  of  six.  As  the 
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pressure  ratio  is  changed  from  six  to  15.2,  the  separation  location,  as  indicated  by  the 
separation  shock,  moves  upstream.  The  trailing  shock  structure,  however,  is  not  well- 
resolved. 

Figures  19  and  20  show  further  the  amount  of  information  provided  by  the  method. 
They  are  plots  of  velocity  vectors.  Separated  areas  are  easily  identified  by  the  abrupt  change 
in  velocity  on  the  afterbody. 

Converged  predictions  of  separation  location  were  achieved  in  approximately  2,100  sec 
of  CRAY-1S  processing  time  for  each  case  presented.  Significant  changes  in  separation 
location  cease  after  about  1,800  sec. 

In  summary,  a  thin-layer,  implicit  Navier-Stokes  method  was  used  to  solve  the  plume- 
induced  separation  problem.  To  achieve  high  static  jet-to-free-stream  pressure  ratios,  the 
nozzle  lip  was  removed  from  the  region  of  the  viscous  solution  and  solved  inviscidly  as  a 
centered  expansion.  The  predictions  were  very  good,  except  at  the  highest  pressure  ratios, 
where  underprediction  of  separation  apparently  resulted  from  poor  grid  resolution.  It  was 
also  demonstrated  that  the  thin-layer  approximation  is  legitimate  for  calculating  large 
separated  How  regions  if  care  is  taken  with  the  turbulence  model,  and  that  the  neai  wake 
could  be  treated  adequately  if  considered  as  being  largely  devoid  of  viscous  effects.  The 
method  demonstrated  a  marked  improvement  in  accuracy  over  the  Chapman-Korst 
component  method,  achieving  a  predictive  accuracy  sufficient  to  allow  the  method  to  be 
used  as  an  engineering  tool.  Also,  it  permits  solutions  of  plume  interference  problems  in  the 
flow  regimes  where  the  Chapman-Korst  component  method  fails.  The  accuracy  of  this 
approach  is  expected  to  improve  with  improved  grid  resolution  when  larger  computer 
memory  becomes  available. 

4.0  CONCLUSIONS  AND  RECOMMENDATIONS 

Two  methods  for  determining  the  extent  of  plume-induced  separation  on  afterbodies  in 
supersonic  flow  have  been  developed.  The  first,  a  Chapman-Korst  component  method  that 
uses  an  external  separation  criterion,  produces  quite  good  predictions  at  very  high  jet-to- 
free-stream  pressure  ratios  where  extensive  separation  is  present.  It  does  not  predict  small 
separation  regions  well,  and  it  fails  when  the  trailing  shock  system  has  embedded  subsonic 
flow.  However,  because  of  its  speed  and  relative  ease  of  use,  coupled  with  the  fact  (hat  it  will 
predict  the  type  of  separation  that  could  blank  a  control  surlace,  it  is  the  first  choice  for 
quick  screening  of  possible  effects  of  chamber  pressure  and  afterbody  design. 

The  second  method  presented  requires  the  solution  of  the  thin-iayer  Navier-Siokcs 
equations.  This  method  shows  its  strength  in  these  very  flow  regimes  where  the  component 
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method  fails,  i.e.,  regimes  where  the  extent  of  separation  is  relatively  small  and  regimes 
where  extensive  subsonic  flow  occurs.  The  thin-layer  form  of  the  equations  proved  to  be 
more  accurate  than  the  component  method  at  low  to  moderate  pressure  ratios.  However,  at 
higher  pressure  ratios,  convergence  becomes  more  difficult  and  careful  attention  to  the  grid 
structure  is  required.  Still,  where  the  component  method  is  a  mature  technology,  the  Navier- 
Stokes  method  shows  potential  for  further  refinement  to  the  point  where  it  may  be  used 
routinely  for  highly  accurate  predictions  of  separated  flows. 

Resources  should  be  turned  now  to  the  development  of  better  Navier-Stokes  methods. 
This  will  require  the  implementation  of  better  turbulence  models,  more  computer  memory, 
and,  of  course,  higher  computational  speeds.  Refinement  of  algorithms  is  also  required  so 
that  the  nozzle-lip  problem  may  be  handled  without  resort  to  artifices  such  as  described  in 
the  text. 
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Figure  2.  Interaction  showing  shear-layer  model. 
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Figure  5.  Afterbody  pressure  distribution,  component  method. 
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Figure  7.  Computational  grid. 
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Figure  8.  Schematic  of  base  region. 
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Figure  9.  Overlapping  grids. 
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Figure  10.  ART  solid-sting  comparison. 
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Figure  12.  Afterbody  pressure  distribution. 
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Figure  13.  Plateau  pressure  comparison. 
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Figure  14.  Base  pressure  comparison 
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APPENDIX  A 

TRANSFORMATION  TO  CURVILINEAR  COORDINATES 

Consider  the  general  form  of  the  dimensionless  conservation  equations  in  Cartesian 
coordinates. 

a,q  +  axE  +  3yF  +  dzG  =  Re"1  (3XR  +  3yS  +  3ZT)  (A-l) 

where  dx  =  3/3x,etc. 

Define  vectors  such  that 


V  =  [q,  E,  F,  G]  =  [V',V2,  V3,  V4] 

W  =  [R,  S,  T]  =  [W,  W2,  W3] 

Equation  (A-l)  may  now  be  written  in  divergence  form  as 

V  -V  =  ~  V  .  W  (A-2) 

Re 

From  tensor  analysis,  viz.  Ref.  23,  in  general  coordinates 

V  .  V  =  J3X,  (J-  1Vi) 

and 

V  •  W  =  Jdxl  ( J  - 1 W*) 

where  the  Jacobian  notation,  J,  is  inverted  to  conform  with  Pulliam  and  Steger.  Elements  V1 
and  W1  are  the  contravariant  components  of  V  and  W.  The  divergence  is  invariant  under 
transformation  as  are  the  vectors;  therefore,  the  only  effort  required  is  in  determining  these 
vector  components  in  the  curvilinear  system. 

Contravariant  vectors  components  transform  from  Cartesian  [t,x,y,z]  to  curvilinear  [r, 
tj,  f]  coordinates  as 


V1  =  3trV> 

V2  =  dt£V*  +  d*£V2  +  3y£V3  +  dz£V4  (A-3) 
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v3  =  a^v1  +  ax>iv2  +  a^v3  +  az)3v4 
v4  =  acfvl  +  axfv2  +  ayrv3  +  az^v4 

In  the  same  way,  W‘  —  W1 
also 

rt  0  0  0 

«l  fx  €y  €* 

J  = 

tfi  Vy  Vz 

ft  fx  fy  fz 

where  £x  =  ax£,  etc. 

From  Eqs.  (A-2)  and  (A-3) 

dr  q  +  a^E  +  dje  +  djG  =  Re-  i(dtR  4-  d„S  +  3rT)  (A-4) 

where 

q  =  J->  V1,  E  =  J-  1  V2,  F  s  J-1  V3,  and  G  =  J”1  V4 

and 

R  ■  J-1  W1,  S  -  J"1  W2,  T  =  J ' 1  W3 

The  form  of  the  conservation  equations  that  is  usually  shown  is  obtained  by  substituting  the 
variables  for  the  vector  components  into  Eq.  (A-4). 

For  example,  consider  the  £  momentum  equation.  The  Cartesian  components  are 

V'  se  q  =  Cu,  V2  =  E  =  eu2  +  p,  V3  =  F  =  qv9,  V4  *  G  =  guw 

W1  =  R  =  tw  W2  =  S  =  T„y,  W3  =  T  == 

These  are  substituted  into  Eqs.  (A-3).  Then  Eqs.  (A-3)  are  substituted  into  Eq.  (A-4). 
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The  left-hand  side  of  Eq.  (A-4)  becomes 

dj-'  e«  aj " 1  £t  eu  aj-1  £x(e«2  +  p)  +  aj~ 1  Sy  (&uv) 

dr  +  d£  &Z 

dS  ~ 1  |z  (guw)  aj~  St  eu 

d£  di} 

+  dJ  ~  1  7?x(eu2  +  p)  aj-'ljyfeuv)  +  d)~'  jz(guw) 
drt  dy  dy 

aj-'Tteu  .  aj-'rx(eu2  +  p)  ,  aj-'iv_(euv)_ 

+  ar  at  '  ar 

+  3J~  'fzCgw)  _  RHS 

dz 

Letting 

U  =  £,  +  £xu  +  iyv  +  £.w 
U  =  T)t  +  rjxu  +  rjyv  +  J?zw 
W  =  ft  +  fxU  +  fyV  +  rzW 

we  get 

dj-1eu  aj~1(QuU  +  £xp)  +  aj~Heuv  +  %p)  +  aj~*(euw  +  fxP)  _  Rj_jg 

+  dt  dy  9t 

and 

RHS  =  Re"1  a£J” '  ($*  7-xx  +  $yT*y  +  ^  +  rjyTxy  +  r/zrxz) 

+  9j-J ~  '(fxTxx  +  i"  yTxy  +  fzTxz) 

or  for  the  thin-layer  approximation: 

RHS  =  Re-1arJ-,(fXTxx  +  tyTjy  +  fz^z) 

The  other  conservation  equations  transform  as  easily. 
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NOMENCLATURE 

a  Sound  speed 

Cp  Static  pressure  coefficient 

cp  Specific  heat  at  constant  pressure 

e  Total  energy 

H  Total  enthalpy 

k  Effective  thermal  conductivity,  k  =  cp  +  Pt/0.9) 

L  Body  length 

£  Length  of  mixing  region;  mixing  length  in  turbulence  model 

4ff  Effective  mixing  length 

fvir  Virtual  mixing  length 

M  Mach  number 

n  Inverse  of  power-law  exponent 

p  Static  pressure 

Pr  Prandtl  numbeT 

R  Radius  from  axis  of  symmetry 

R  Average  radius 

Reu  Unit  Reynolds  number 

U  Velocity  of  high-speed  edge  of  shear  layer  in  direction  of  shear  layer 

u  Velocity  component  in  Cartesian  x  direction;  velocity  in  direction  of 

development  shear  layer 
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U,  V,  W  Contravariant  velocity  components 

v,  w  Velocity  component  in  Cartesian  y  and  z  directions,  respectively 


X 

x 


y 

e 


5 

y 


% 


<*> 

T 


fit 


Coordinate  along  slipline  in  recompression  region 

Cartesian  axial  coordinate  or  mixing  coordinate  in  direction  of  developing 
mixing  layer 

Cartesian  coordinate  transverse  to  x  coordinate 

Angle  between  slipline  and  mixing  layer,  Fig.  3;  also  used  for  nozzle  and 
afterbody  angles 

Boundary-layer  thickness 

Ratio  of  specific  heats 

Curvilinear  coordinate  approximately  normal  to  body,  also  boundary-layer 
coordinate,  y/5 

Dimensionless  mixing  variable;  body-conforming  curvilinear  coordinate  in 
azimuthal  direction 

Dimensionless  position  parameter,  ab/t 
Velocity  ratio  in  shear  layer,  u/U 
Shear  stress;  time 

Effective  viscosity:  molecular  +  turbulent,  except  when  subscripted 
Turbulent  viscosity 


Molecular  viscosity 

Body-conforming  curvilinear  coordinate  in  flow  direction 


Density 
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a  Mixing  or  spreading  parameter 

a0  Incompressible  a,  a  reference  value 

w  Vorticity 

SUBSCRIPTS 

1  Beginning  of  recompression;  also  upstream  of  shock  wave 

2  End  of  recompression;  also  downstream  of  separating  shock  wave 

B  Separated  region;  base  region 

D  Dividing  streamline 

E  Nozzle  exit  condition 

oo  Free-stream  condition 

L  Low-speed  edge  of  mixing  layer 

M  Molecular 

m  Inviscid  boundary  with  respect  to  origin  of  mixing  profile 

N  Nozzle 

r  Downstream  of  recompression  shock  wave 

S  Stagnating  streamline 

51  Stream  1  (jet  plume) 

52  Stream  2  (outer  stream) 

T  Total  condition 

U  High-speed  edge  of  mixing  layer 
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